#PBS -N CHG002591
#PBS -l nodes=1:ppn=20
#PBS -l mem=120G
#PBS -q dna
cd $PBS_O_WORKDIR


# mkdir -p .//2.mapping;
#bwa mapping ----
# /mnt/ilustre/app/rna/bin/bwa aln -e 3 -i 10 -t 10 -R 100 -q 20 /mnt/ilustre/app/medical/.data/ref/b37/b37.fa ../hg19/1.QC/CHG002591_trim1.fastq > .//2.mapping/CHG002591-1.sai;
# /mnt/ilustre/app/rna/bin/bwa aln -e 3 -i 10 -t 10 -R 100 -q 20 /mnt/ilustre/app/medical/.data/ref/b37/b37.fa ../hg19/1.QC/CHG002591_trim2.fastq > .//2.mapping/CHG002591-2.sai;
# /mnt/ilustre/app/rna/bin/bwa aln -e 3 -i 10 -t 10 -R 100 -q 20 /mnt/ilustre/app/medical/.data/ref/b37/b37.fa ../hg19/1.QC/CHG002591_trim_unpaired.fastq > .//2.mapping/CHG002591-s.sai;

# /mnt/ilustre/app/rna/bin/bwa sampe -a 500 -r "@RG\tID:LS5993870\tSM:CHG002591\tLB:CHG002591\tPU:run barcode\tPL:Illumina\tDS:resequencing\tCN:MajorBio" /mnt/ilustre/app/medical/.data/ref/b37/b37.fa .//2.mapping/CHG002591-1.sai .//2.mapping/CHG002591-2.sai ../hg19/1.QC/CHG002591_trim1.fastq ../hg19/1.QC/CHG002591_trim2.fastq | /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools view -bS - > .//2.mapping/CHG002591-sampe.bam;
# /mnt/ilustre/app/rna/bin/bwa samse -r "@RG\tID:LS5993870\tSM:CHG002591\tLB:CHG002591\tPU:run barcode\tPL:Illumina\tDS:resequencing\tCN:MajorBio" /mnt/ilustre/app/medical/.data/ref/b37/b37.fa .//2.mapping/CHG002591-s.sai ../hg19/1.QC/CHG002591_trim_unpaired.fastq | /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools view -bS - > .//2.mapping/CHG002591-samse.bam;

# /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools sort -m 3000000000 .//2.mapping/CHG002591-sampe.bam .//2.mapping/CHG002591-sampe.sort;
# /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools sort -m 3000000000 .//2.mapping/CHG002591-samse.bam .//2.mapping/CHG002591-samse.sort;
# /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools merge .//2.mapping/CHG002591.all.sort.bam  .//2.mapping/CHG002591-sampe.sort.bam .//2.mapping/CHG002591-samse.sort.bam;
# /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools index .//2.mapping/CHG002591.all.sort.bam;

:<<!
#picard rmdup ----
java -Xmx40G -jar /mnt/ilustre/users/jie.chen/software/picard-tools-1.102/MarkDuplicates.jar TMP_DIR=.//2.mapping/tmp MAX_FILE_HANDLES=1000 VALIDATION_STRINGENCY=SILENT ASSUME_SORTED=true REMOVE_DUPLICATES=true I=.//2.mapping/CHG002591.all.sort.bam O=.//2.mapping/CHG002591.all.sort.rmdup.bam M=.//2.mapping/CHG002591.metrics 2> .//2.mapping/CHG002591.rmdup.log;
/mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools index .//2.mapping/CHG002591.all.sort.rmdup.bam;
!

#GATK realign ----
java -Xmx30G -jar /mnt/ilustre/users/jie.chen/tools/GenomeAnalysisTK.jar -T RealignerTargetCreator -R /mnt/ilustre/app/medical/.data/ref/b37/b37.fa -I .//2.mapping/CHG002591.all.sort.rmdup.bam -log .//2.mapping/CHG002591.intervals.log -o .//2.mapping/CHG002591.realign.intervals;

java -Xmx30G -jar /mnt/ilustre/users/jie.chen/tools/GenomeAnalysisTK.jar -T IndelRealigner -R /mnt/ilustre/app/medical/.data/ref/b37/b37.fa -I .//2.mapping/CHG002591.all.sort.rmdup.bam -log .//2.mapping/CHG002591.realign.log -targetIntervals .//2.mapping/CHG002591.realign.intervals -o .//2.mapping/CHG002591.all.sort.rmdup.realign.bam;


#samtools mapping stats ----
/mnt/ilustre/users/jie.chen/software/circos_script/iTools Fatools stat -InPut /mnt/ilustre/app/medical/.data/ref/b37/b37.fa -OutPut .//2.mapping/chr.list;
perl /mnt/ilustre/users/jie.chen/software/circos_script/Mapping_Info.pl /mnt/ilustre/users/jie.chen/software/samtools-0.1.18/samtools .//2.mapping/CHG002591.all.sort.rmdup.bam .//2.mapping/chr.list .//2.mapping/CHG002591.mapstat.tvs;


# mkdir -p .//3.varcalling;
#GATK snp and indel calling ----
java -Xmx30G -jar /mnt/ilustre/users/jie.chen/tools/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /mnt/ilustre/app/medical/.data/ref/b37/b37.fa -I .//2.mapping/CHG002591.all.sort.rmdup.realign.bam -o .//3.varcalling/CHG002591.snp.vcf -stand_call_conf 50 -stand_emit_conf 50 -dcov 200 -glm SNP

java -Xmx30G -jar /mnt/ilustre/users/jie.chen/tools/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /mnt/ilustre/app/medical/.data/ref/b37/b37.fa -I .//2.mapping/CHG002591.all.sort.rmdup.realign.bam -o .//3.varcalling/CHG002591.indel.vcf -stand_call_conf 50 -stand_emit_conf 50 -dcov 200 -glm INDEL


#vcfutils varFilter ----
/mnt/ilustre/app/dna/other_software/linaiting/tools/Module_Reseq/samtools-0.1.18/vcfutils.pl varFilter -d 10 -D 10000 .//3.varcalling/CHG002591.snp.vcf > .//3.varcalling/CHG002591.snp.filt.vcf;
/mnt/ilustre/app/dna/other_software/linaiting/tools/Module_Reseq/samtools-0.1.18/vcfutils.pl varFilter -d 10 -D 10000 .//3.varcalling/CHG002591.indel.vcf > .//3.varcalling/CHG002591.indel.filt.vcf;


# rm .//2.mapping/CHG002591-1.sai .//2.mapping/CHG002591-2.sai .//2.mapping/CHG002591-s.sai;
# rm .//2.mapping/CHG002591-sampe.bam .//2.mapping/CHG002591-sampe.sort.bam .//2.mapping/CHG002591-samse.bam .//2.mapping/CHG002591-samse.sort.bam;
